This file explains how we ran the statisticla tests on the aperture size measurements. It is based on the Aps dataframe, which we created in the Basic Shell Information.R file. It is worth to redo these steps because we won’t summarise it as much this time. First, we need to source the ‘Farasan Database.R’ file through:
source('Farasan Database.R')
Then we extract the aperture measurements from the Apertures list:
Aps <- data.frame()
f <- data.frame()
for (i in 1:length(Apertures)) {
for (j in 1:length(Apertures[[i]])) {
if ("Strombus.fasciatus..............Born.1778..Aperture" %in% names(Apertures[[i]][[j]]) == FALSE) {
next
}
f <-
data.frame(c(as.character(names(Apertures[i]))), c(as.character(names(Apertures[[i]][j]))), Apertures[[i]][[j]][["Strombus.fasciatus..............Born.1778..Aperture"]])
colnames(f) <- c('Site', 'Sample', 'Aperture')
Aps <- rbind(Aps, f)
}
}
rm(f, i, j, Areas, Sites, Species)
To get a simple list of each measurement grouped by site and bulk sample number:
Aps[1:3,1:3]
## Site Sample Aperture
## 1 KM1048 926 24.02
## 2 KM1048 926 18.13
## 3 KM1048 926 20.07
To look at the general statistics for each areas, we added a column to Aps containing the respective bay areas and summarise the dataframe using the summarySE function:
Aps$Area <- c(rep("Khur Maadi",4122), rep("Janaba West",6171),rep("Janaba East",5169))
source("summarySE.R")
Aps_area <-summarySE(na.omit(Aps), measurevar = "Aperture", groupvars = c("Area"))
## Area N Aperture_mean Aperture_median sd
## 1 Janaba East 5169 19.67116 19.60 1.483827
## 2 Janaba West 6171 21.39223 21.30 1.648272
## 3 Khur Maadi 4122 22.87554 22.85 1.821071
Doing an one-way analysis of variance (anova) will show whether there are statistically significant differences between the means of each bay area.
Aps_area_aov <- aov(Aperture ~ Area, data = Aps)
summary(Aps_area_aov)
## Df Sum Sq Mean Sq F value Pr(>F)
## Area 2 23880 11940 4415 <2e-16 ***
## Residuals 15459 41808 3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We then applied a post-hoc tukey test, to understand between which areas the statistically significant differences lie.
TukeyHSD(Aps_area_aov)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Aperture ~ Area, data = Aps)
##
## $Area
## diff lwr upr p adj
## Janaba West-Janaba East 1.721064 1.648385 1.793743 0
## Khur Maadi-Janaba East 3.204381 3.123888 3.284873 0
## Khur Maadi-Janaba West 1.483317 1.405778 1.560856 0
To illustrate these differences, our paper used raincloud plots, as well as other design-related packages, such as viridis and cowplot.
library(ggplot2)
library(cowplot)
library(viridis)
source('R_rainclouds.R')
Fig6a <- ggplot(Aps, aes(
x = Area,
y = Aperture,
fill = Area,
colour = Area
)) +
geom_flat_violin(
position = position_nudge(x = +0.155, y = 0),
adjust = 0.75,
trim = FALSE,
alpha = .99
) +
geom_point(position = position_jitter(width = .1),
size = 0.05,alpha=0.2
) +
geom_boxplot(
position = position_nudge(x = +0.155, y = 0),
aes(x = Area, y = Aperture,colour=Area),
outlier.shape = NA,
alpha = 0.3,
width = .1,
colour = "BLACK",show.legend = FALSE
) +
ylab('Aperture size') + xlab('') + coord_flip() + guides(fill = FALSE, colour = FALSE) +
scale_color_manual(values=c(viridis_pal()(60)[10],viridis_pal()(6)[4],viridis_pal()(60)[59]))+
scale_fill_manual(values=c(viridis_pal()(60)[10],viridis_pal()(6)[4],viridis_pal()(60)[59]))+
theme_cowplot()+
theme(axis.text.y.left = element_text(vjust=-2))+theme(axis.ticks.y=element_blank())
Fig6a
To do the same tests between sites, we have to take similar steps, only this time we are also using the Site column to group aperture measurements.
##Kruskal Wallis The same data on site level reveals that there are variations in sample size, normal distributions and variance between the sites, which are all assumptions for the ANOVA test used above. We will thus switch to the non-parametric Wilcoxon rank-sum test and use the pairwise.wilcox.test() function
PW_Aps <- pairwise.wilcox.test(Aps$Aperture, Aps$Site, p.adjust.method = "bonferroni")
PW_Aps
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: Aps$Aperture and Aps$Site
##
## KM1048 KM1050 KM1051 KM1052 KM1053 KM1054 KM1056 KM1057
## KM1050 1.00000 - - - - - - -
## KM1051 1.00000 0.60812 - - - - - -
## KM1052 1.00000 0.02953 1.00000 - - - - -
## KM1053 1.00000 1.00000 1.00000 1.00000 - - - -
## KM1054 1.00000 1.00000 1.00000 0.60836 1.00000 - - -
## KM1056 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 - -
## KM1057 1.00000 1.00000 1.00000 0.15740 1.00000 1.00000 1.00000 -
## KM1304 1.00000 0.01776 1.00000 1.00000 1.00000 1.00000 1.00000 0.10905
## KM1307 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## KM1313 1.00000 1.00000 0.21637 0.00616 1.00000 1.00000 1.00000 1.00000
## KM1317 1.00000 1.00000 1.00000 0.67309 1.00000 1.00000 1.00000 1.00000
## KM1328 1.00000 1.00000 1.00000 0.19374 1.00000 1.00000 1.00000 1.00000
## KM1330 1.00000 1.00000 1.00000 0.11719 1.00000 1.00000 1.00000 1.00000
## KM1335 1.00000 1.00000 1.00000 0.20132 1.00000 1.00000 1.00000 1.00000
## KM1336 1.00000 1.00000 1.00000 0.48361 1.00000 1.00000 1.00000 1.00000
## JW1705 1.00000 7.7e-06 1.00000 1.00000 1.00000 0.52703 1.00000 0.00011
## JW1727_A 1.00000 < 2e-16 0.00331 1.00000 3.7e-16 0.00090 1.00000 < 2e-16
## JW1727_B 1.00000 < 2e-16 0.00051 1.00000 < 2e-16 0.00047 1.00000 < 2e-16
## JW1807 1.00000 1.2e-09 1.00000 1.00000 1.00000 1.00000 1.00000 < 2e-16
## JW1864 1.00000 < 2e-16 1.5e-05 0.57737 < 2e-16 0.00013 1.00000 < 2e-16
## JW2298 1.00000 < 2e-16 0.04956 1.00000 9.3e-12 0.00231 1.00000 < 2e-16
## JW3120 1.00000 < 2e-16 1.00000 1.00000 2.0e-07 0.02067 1.00000 < 2e-16
## JE0003 1.00000 0.00153 0.19107 1.00000 0.00407 0.19211 1.00000 0.00202
## JE0004 1.00000 < 2e-16 < 2e-16 1.2e-15 < 2e-16 3.4e-10 1.00000 < 2e-16
## JE0078 1.00000 < 2e-16 < 2e-16 4.1e-10 < 2e-16 2.5e-08 1.00000 < 2e-16
## JE0087 1.00000 < 2e-16 5.0e-16 1.1e-06 < 2e-16 3.9e-07 1.00000 < 2e-16
## JE5642 1.00000 < 2e-16 5.3e-08 0.00568 < 2e-16 2.1e-05 1.00000 < 2e-16
## KM1304 KM1307 KM1313 KM1317 KM1328 KM1330 KM1335 KM1336
## KM1050 - - - - - - - -
## KM1051 - - - - - - - -
## KM1052 - - - - - - - -
## KM1053 - - - - - - - -
## KM1054 - - - - - - - -
## KM1056 - - - - - - - -
## KM1057 - - - - - - - -
## KM1304 - - - - - - - -
## KM1307 1.00000 - - - - - - -
## KM1313 0.00520 1.00000 - - - - - -
## KM1317 0.58399 1.00000 1.00000 - - - - -
## KM1328 0.11911 1.00000 1.00000 1.00000 - - - -
## KM1330 0.10090 1.00000 1.00000 1.00000 1.00000 - - -
## KM1335 0.20101 1.00000 1.00000 1.00000 1.00000 1.00000 - -
## KM1336 0.47794 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -
## JW1705 1.00000 1.00000 6.4e-06 1.00000 0.01718 0.06750 0.13757 0.67647
## JW1727_A 1.00000 0.00028 < 2e-16 1.1e-10 8.1e-12 2.8e-10 2.4e-13 1.9e-11
## JW1727_B 1.00000 9.5e-05 < 2e-16 1.4e-11 3.0e-12 7.0e-11 7.1e-14 2.3e-12
## JW1807 1.00000 1.00000 8.2e-08 1.00000 0.01638 0.10223 0.07202 0.84577
## JW1864 1.00000 5.6e-06 < 2e-16 1.9e-14 5.5e-15 4.5e-13 < 2e-16 2.4e-15
## JW2298 1.00000 0.00087 < 2e-16 5.0e-09 2.7e-09 4.1e-09 1.2e-10 1.4e-09
## JW3120 1.00000 0.04652 < 2e-16 7.9e-06 3.0e-07 1.8e-06 1.1e-07 2.5e-06
## JE0003 1.00000 0.04016 0.00086 0.00382 0.02436 0.00269 0.00533 0.00465
## JE0004 3.6e-16 3.3e-15 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16
## JE0078 9.2e-10 1.2e-12 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16
## JE0087 5.3e-06 6.1e-11 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16
## JE5642 0.01192 4.5e-07 < 2e-16 6.5e-15 4.1e-15 3.6e-13 < 2e-16 1.0e-15
## JW1705 JW1727_A JW1727_B JW1807 JW1864 JW2298 JW3120 JE0003
## KM1050 - - - - - - - -
## KM1051 - - - - - - - -
## KM1052 - - - - - - - -
## KM1053 - - - - - - - -
## KM1054 - - - - - - - -
## KM1056 - - - - - - - -
## KM1057 - - - - - - - -
## KM1304 - - - - - - - -
## KM1307 - - - - - - - -
## KM1313 - - - - - - - -
## KM1317 - - - - - - - -
## KM1328 - - - - - - - -
## KM1330 - - - - - - - -
## KM1335 - - - - - - - -
## KM1336 - - - - - - - -
## JW1705 - - - - - - - -
## JW1727_A < 2e-16 - - - - - - -
## JW1727_B < 2e-16 1.00000 - - - - - -
## JW1807 1.00000 < 2e-16 < 2e-16 - - - - -
## JW1864 < 2e-16 0.00316 1.00000 < 2e-16 - - - -
## JW2298 3.7e-10 1.00000 1.00000 < 2e-16 0.35921 - - -
## JW3120 9.5e-05 1.4e-05 1.9e-07 2.0e-12 < 2e-16 0.58862 - -
## JE0003 0.01027 1.00000 1.00000 0.01014 1.00000 1.00000 0.38053 -
## JE0004 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 0.63341
## JE0078 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.00000
## JE0087 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 < 2e-16 1.00000
## JE5642 < 2e-16 8.3e-08 0.00032 < 2e-16 0.00197 3.8e-06 1.5e-15 1.00000
## JE0004 JE0078 JE0087
## KM1050 - - -
## KM1051 - - -
## KM1052 - - -
## KM1053 - - -
## KM1054 - - -
## KM1056 - - -
## KM1057 - - -
## KM1304 - - -
## KM1307 - - -
## KM1313 - - -
## KM1317 - - -
## KM1328 - - -
## KM1330 - - -
## KM1335 - - -
## KM1336 - - -
## JW1705 - - -
## JW1727_A - - -
## JW1727_B - - -
## JW1807 - - -
## JW1864 - - -
## JW2298 - - -
## JW3120 - - -
## JE0003 - - -
## JE0004 - - -
## JE0078 < 2e-16 - -
## JE0087 < 2e-16 < 2e-16 -
## JE5642 < 2e-16 4.5e-08 1.00000
##
## P value adjustment method: bonferroni
Similar to the graph above, we illustrate the site differences using raincloud plots, as well as other design-related packages, such as viridis and cowplot. This time, the sites are coloured by bay area, however they are in the order of largst shells (top) to smallest shells (bottom). The general trend, however, remains and sites usually stick to their bay area.
library(ggplot2)
library(cowplot)
library(viridis)
source('R_rainclouds.R')
Fig6b <- ggplot(Aps, aes(
x = reorder(Site, Aperture, FUN = mean,order=TRUE),
y = Aperture,
fill = Area,
colour = Area, order=Aperture
)) +
geom_flat_violin(
position = position_nudge(x = +0.2, y = 0),
adjust = 0.75,
trim = FALSE,
alpha = .99
) +
geom_point(position = position_jitter(width = 0.1),
size = 0.5,alpha=0.3
) +
geom_boxplot(
position = position_nudge(x = +0.2, y = 0),
aes(x = Site, y = Aperture,colour=Area),
outlier.shape = NA,
alpha = 0.3,
width = .1,
colour = "BLACK"
) +
ylab('Aperture size') + xlab("") + coord_flip() +guides(fill = guide_legend(reverse = TRUE),color = guide_legend(reverse = TRUE))+
scale_color_manual(values=c(viridis_pal()(3)[1],viridis_pal()(6)[4],viridis_pal()(60)[59]))+
scale_fill_manual(values=c(viridis_pal()(3)[1],viridis_pal()(6)[4],viridis_pal()(60)[59]))+ theme_cowplot()+ ggtitle("")+
theme(axis.text.y.left = element_text(vjust=-0))+
theme(axis.ticks.y=element_blank())+
theme(legend.position="right")
Fig6b
Finally, we want to integrate both raincloud plots (Areas and Sites) into one figure. For this we use the patchwork package.
library(patchwork)
Fig6 <- Fig6a+Fig6b
Fig6 + plot_annotation(
title = 'Aperture sizes',
subtitle = 'Differences between bay area on the left and differences between sites on the right')
setwd("~/Farasan-Shell-Data-master/Figures and supplementary material")#Remember to change the file path
ggsave("Figure_6.png", height = 200, width = 300, units = "mm")